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ABSTRACT 

We analyze the distribution of stars of arbitrary mass function ^(m) around a massive black hole 
(MBH) . Unless ^ is strongly dominated by light stars, the steady-state distribution function approaches 
a power-law in specific energy x = —E/ma'^ < Xmax with index p — m/AMo, where E is the energy, 
a is the typical velocity dispersion of unbound stars, and Mq is the mass averaged over m^xP^^^. For 
light-dominated ^, p can grow as large as 3/2 - much steeper than previously thought. A simple 
prescription for the stellar density profile around MBHs is provided. We illustrate our results by 
applying them to stars around the MBH in the Milky Way. 

Subject headings: Galaxy: kinematics and dynamics — stellar dynamics — black hole physics 



1. INTRODUCTION 

A massive black hole (MBH) of mass M, dominates 
the dynamics of stars within its radius of influence = 
GM,/a^ = 2.3(A/./3 x 106Mq)(ct/75 km s~^)-^ pc, 
where a is the typical star velocity dispersion at r > 
and values quote d correspond to the Milky Way center 
(| Alexander! 120051 ). The steady-state distribution func- 
tion (DF) of such stars, first derived (for simple, dis- 
crete stellar mass functions, MFs) by Bahcall and Wolf 
(1976, 1977; henceforth BW76,77) is useful in the study 
of galactic centers and possibly globular clusters. In 
the simple case where all stars have equal mass, BW76 
showed that an n cx r~'^/'^ density cusp forms; a some- 
what flatter cusp was subsequently ide ntified around 
SgA* (|Alexandeilll999l : IGenzel et al.li2003l : ISchodel et all 
[20071) . 

Significant mass segregation is expect ed for realis- 
tic MFs (iMirald a-Escudc and Gould 2000i: ISchodel et all 
120071: lO'Leary et all 120081 and references therein), as 
dynamical friction slows heavy stars that sink to- 
wards the MBH while light stars are pushed out- 
wards. This effect is essential in modeling var- 
ious physical processes, s u ch as tidal disruption 
jLightman and Shapiro' '1973 iMagorrian and Tremain^ 
[1999: Syer and Ulmer 199 9) and grav it ational wave 
emiss i on ( Hopman and Ale xander! 120061 : iFreitag et al.l 
[20061 : [Hopman et al., i2007.) . and may be used for ex- 
ample to test if intermedi ate-mass b lack holes ex- 
ist in globular clusters CGill et_alj [2008[ ). Re- 
cently, lAlixander^nTJoEman (200!jlienceforth AH09) 
showed strong mass segregation in the limit of a strongly 
bimodal MF when the mass ratio is large and massive 
stars are rare. 

In spite of extensive numerical studies of mass segrega- 
tion around a MBH (AH09 and references therein), an- 
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alytical modeling (BW77) is limited to simple, discrete, 
and extreme MFs where the most massive species consti- 
tutes more than half of the star population. Most results 
are limited to few stellar species and to a narrow mass 
range. In this letter we analytically study the steady- 
state single-star DF / within r < rh, for an arbitrary, 
continuous MF. We confirm our results numerically for 
a wide range of MFs. 

Following BW76,77, we assume (i) spherical spatial 
symmetry; (ii) isotropic velocities; (iii) Keplerian orbits; 
(iv) binaries are negligible; (v) small angle, uncorrelated, 
local scatterings; (vi) loss cone effects can be neglected; 
(vii) isothermal distribution of unbound stars with tem- 
perature /icr^; (viii) stars are destroyed when their spe- 
cific energy drops below a threshold x = —E/ma^ = 
Xmax, with E the energy. 

2. EQUAL MASS STARS 

Consider the case where all stars have mass m. The 
dimensionless diffusion rate of stars through x is (BW76) 
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The boundary conditions are /(x < 0) = e^'"/^ and 
fix > Xmax) — 0. Eq. ll]) with these boundary condi- 
tions uniquely determines / and Q as shown below. The 
spatial number density of stars is given by (BW76) 
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Eq. (HI) describes the balance between diffusion towards 
the MBH and replenishment by new stars. BW76 found 
that Q is essentially determined by the diffusion rate at 
a bottleneck near Xmax, where the replenishment rate 
diminishes. As Q — ^/xmax is typically small, one can 
approximately solve Eq. ^ by setting (3 = 0. 

We convert Eq. ^ to an ordinary differential equation 
(ODE) by repeated operations of differentiation with re- 
spect to x and isolation of integral terms. This yields a 
nonlinear, fourth-order ODE for /, 
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Here we defined the local power-law index of /, p 
din //din a; = xw, the local power-law index of w 
(In/)', q = dlnw/dlnx, and the operators 

A{f, x) = V - 5g - 2x'^w"/w , 
B{f, x) = 5q- 12x^w"^/{w'w) + Sx^w^^^ /w 



(4) 
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where D{f,x) = q+p— 5/2. The boundary conditions at 
a; = fix / and /' there, so the non-exponential solution 
f{x) is completely determined for given Q. A unique 
value of Q guarantees that / vanishes (for the first time) 
at Xmax, proving the uniqueness of the steady-state so- 
lution (an "exercise left for the reader" by BW76). 
For a power-law DF / = fox^" , Eq. ^ becomes 

3/2 {po - 1/4) = Qx'/^-'P^po {po - 1) (Po - 5/2) . (5) 

This compactly reproduces the BW76 results: although 
a finite, ener gy-independent flow requires po — 3/4 
(jPeeblesI I1972L implying a flow of stars away from the 
MBH, see BW76), a steady-state with po = 1/4 is closer 
to the actual DF because Q is negligibly small. Note that 
the power-law assumption becomes inconsistent at high 
energies, where Qx^^^^^p° oc a; is large. 

Henceforth we assume Q — 0, justified for single mass 
stars if Xmax S> 1. An exponential DF with q — solves 
Eq. fSj, but does not satisfy the boundary conditions. 
For q 0, Eq. ([3]) becomes a third order nonlinear ODE, 



B + 4pA = 
with general solution of the form 
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+ C3 + CiCsWix; Cl, C2, ca) , (7) 



with ci,C2,C3 constants. 

Useful results can be deduced directly from Eqs. ^ 
and ([7|), without determining the analytic properties of 
W. Whenp is much smaller (larger) than 1/4, Eq. ^ be- 
comes approximately i? ~ (A ~ 0), which corresponds 
to taking W (taking C3) to zero in Eq. ([7]). In this case 
f{x ;» 1) oc exp[J^ w{x')dx'] is either nearly constant, 
or (super-) exponential in x. In order to maintain rea- 
sonable stellar densities, there must be an approximate 
balance between the two terms, implying that p ~ 1/4 
far from the boundaries. Note that for p — 1/4, the six 
terms of A and B in Eq. ([4|) precisely cancel in pairs. 

For small x, where < p ^ xw ^ 1, Eq. ^ becomes 
B ~ 0; this corresponds to Eq. ([7]) with W 0, and 
C3 = m//i according to the boundary conditions at a; = 0. 
At intermediate energies where approximately / oc xP°, 
Eq. ^ becomes -12w"'^ + Sw'^'w' ~ 0, with all other 
terms smaller by factors > 5 (> 10) for \pq\ < 1/2 {\pq\ < 
1/4), so u; ~ po{x + xi)^^ + X2^ . At high energies near 
the disruption energy, p is large and negative, Eq. ^ 
becomes A ~ and so C3 — 0; the boundary condition 

at Xmax then yields C2 = c\/ (3n) = Xmdx with n S N. 
Smoothly combining these approximate solutions yields 



an approximation for the equal mass case, 
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if X > Xta, 
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with constants X2 , Xta and / determined for example by 



a;.. 



continuity of /, SnXm^Jx / {x - x^Jx) = Pqx ^ 
and its first derivative evaluated at the turn-around en- 
ergy Xta- The results depend weakly on n; it is nat- 
ural to fix n — 1. The previous paragraph suggests 
that Po — 1/4 as Xmax increases. Numerically solv- 
ing Eq. (P) yields po = {'^/^)p[xmax), where the cor- 
rection factor p ~ 1 approaches unity as Xmax becomes 
large, e.g., p(lO^) = 1.18 and p(lO^) = 1.08. It also 
depends weakly on the energy range used to measure p 

(here xUmx < x < xHax, henceforth). Figured] illustrates 
the approximation. 




and /x/m = 1, according to (i) exact (Q ~ 8 X 10 *, solid) and (ii) 
approximate (Q = 0, dashed) numerical solutions of Eq. |[TJ, (iii) 
Eq. ((Sj (dash-dotted), and (iv) / ^ 2x^^'^ (dotted). Logarithms 
are base 10 and convergence is better than 1%, henceforth. 



3. CONTINUOUS MASS FUNCTION 

Consider stars with variable mass m in some range 
Ml < to < Mh. We generalize the discrete, multiple- 
mass version of Eq. ^ (BW77) to the continuum limit. 



Mt, 



dx' [Max(a;, a;')] 
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Q{m)= / m' dm' 

J Ml 

X [m/(a;, m)dx'f[x', m') - m'dxf[x, m)f{x', to')] , 

where / is now the DF in x — m phase space. Assuming 
unbound stars with MF i^(to), f{x < 0, m) = ^{m)e'^"^^'^ . 
The density nm(r) is related to f{x, m) as in Eq. ([2]), such 
that nm{r > rt) ~ i{m). 

3.1. Negligible flow 

Consider the negligible flow limit Q ^ 0, which holds if 
the MF is not strongly dominated by light stars (AH09), 
as shown in ii3.2l Dividing Eq. ([9]) by mf{x, m) and dif- 
ferentiating with respect to m we find dmii'mD^^dxf] = 
0, implying that / must have the functional form 

/(x, to) = ^{m)h{x)"' = ,e("^)e"''^''^ , (10) 

with k = \nh. Eq. (fTO)) implies that the local power- 
law index of /, p{x,m) = dhi f /dlnx = mxk'{x), is 
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linear in m for fixed x. Tliis result, valid for negligible Q, 
generalizes the BW77 result p{x, mi) /mi = p{x, m2)/m2 
to a continuous MF. 

Setting Q = in Eq. IH), using Eq. (fTO|) . and repeat- 
edly differentiating with respect to x and isolating inte- 
gral terms, yields an ODE which, for q 7^ 0, becomes 



B + {p/P)A ^ 



(11) 



Here, p, A and B are functionals of /(m, x) and x, defined 
as in the equal mass case {e.g., Eq. ([4])), and we defined 



P{x, m) \ 
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Ml < A/o = (™") //(m) / = (m),^/ < Mh , (13) 

with = / X(m)Y{m)dm/ jY{m)dm averaging. 

The full distribution may be found by solving Eq. (|lip 
for k{x) under the boundary conditions k{x < 0) = x/^ 
and k ^ 00 as x ^ •^max- 

Note that P{x,Ml) < 1/4 and P{x,Mh) > 1/4. If 
the mass range is sufficiently narrow or ^ is sufficiently 
dominated by high masses such that P{x, Mh) — 1/4, we 
recover the equal mass case Eq. ^ for the most massive 
stars. The full distribution then becomes 



/(x,m;M)~f(m)/e,(a:;M)™/*'" , 



(14) 



with /eg(a;;/i) the DP of equal-mass stars (Figure[T]). 

In the general case we may proceed as in [J2j / will 
be constant, vanish or diverge unless the two terms in 
Eq. pT|) are approximately balanced, i.e., p ~ P. If P 
depends only weakly on x, this implies a power-law DP, 
f (X x^ . Thus, stars with mass equal to the weight- 
averaged mass Mq tend to have a p = 1/4 cusp as in the 
equal mass case. Taking some typical energy x^^^ with 
a < 1 in Eq. (|13p . we may calibrate a against numerical 
solutions of Eq. This yields a ~ 1 for a wide range 
of MPs tested. Hence 



(15) 



(16) 



p{m) ~ p{x 

max )m/4Afo 
and the full distribution is approximately 

fix, to; fi) ~ ^{m)feg {x; p)^'^'^" 

Eqs. (jlSp and agree (to better than a factor of 2 in 
/) with numerical solutions of Eq. ^ for various MPs, 
such as the discrete MPs tested by BW77 and AH09, the 
Salpeter MF (see Figure [2|) and a wide range of power- 
law MPs (see Figure O. 
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Fig. 2.— Unevolved Salpeter MF § oc Af-2-35 ^j^h Mh/Mj^ = 

100, ^l = {M^MhY^'^ and Xmax = 10*. DF shown (top to bot- 
tom) for m/Mij = 100, 67, 34 and 1.5, found by numerically solving 
Eq. lO (solid) and from Eq. III6I I (dashed). Also shown for refer- 
ence are power law curves 

6a;i/2 and 2x^1^^ (dotted). 



Eq. (fT5)) implies that for any MF, p{m) is a decreasing 
function of Xmax- In particular, the power-law index of 
the most massive species, pn, asymptotically approaches 
1/4 as Xmax ^ 00, so / approaches Eq. (fT4| . In cases 
where we may approximate / '--^ ^ in Eq. (|13p . for ex- 
ample if TO^ is strongly peaked in a mass range where 
j l^ — h™ varies little, we have Mq ~ {m)m^. 

For concreteness, consider unbound stars with a power- 
law MP ^(to) c>c to" in some mass range Ml < to < 
Mh — CMl- In this case we recover Eq. ([TT|) with 



P = - 



Mk r(2 + a, -Mnk, -MLk) 
~ r{3 + a,-MHk,-MLk) 
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and r(a,6,c) = J^; f"- 

The power-law index p{m) may be found from Eq, 
^lns-"P,^ln 



* dt the incomplete F-function. 
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(for one mass, then using p cx to). Indices pl, Ph cal- 
culated from Eq. are illustrated in Figure [S] in the 
^ — a plane. For small (very negative) a we may estimate 
p using the fc — )■ limit where Mq ^ (to)„j, whence ^ 
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Fig. 3.— Power-law MFs ^{Ml < m < Mh = (Ml) oc m°' 
in the — a phase space for Xmax = 10* and = (MlMjj)^^^ ■ 
Shown are spectral indices Pl,Ph of the lightest, heaviest (color 
scale, black contours) stars, found by numerically solving Eq. 
The power-law assumption f{x, Mh) oc x^i^ fails (p-value of fit 
larger than 1/2, henceforth) in the region enclosed by red contours. 
The linear scaling p oc m breaks-down beneath the dashed yellow 
contour (wavy appearance of these contours is a resolution effect). 
Approximation Eq. I I15I I f or p n is shown (dashed blue contours) 
outside the <5 < (see Eq. | |23| |) region, which is enclosed by thick 
blue contours. 



Special cases: P£f (a = —2) : 



4{C-1) 



and p_ff(Q! = —3) 



C-i 

41nC' 
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3.2. Non-Negligible flow 

When Q cannot be neglected, we may eliminate the x' 
integral in Eq. ([9]) by repeated operations of differentia- 
tion with respect to x and isolation of integral terms. If 
we assume, in addition, a power-law energy dependence 



/ 



_-Mh 



, we arrive at a generalization of Eq. (O, 
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with abbrev. p = p{m) and p' = p{m'). For Q = Q and 
p (X m, this reproduces the results of p = m/4Mo. 

Eq. (pOj) suggests that p < 3/2 in all cases, other- 
wise the left hand side becomes large and strongly x- 
dependent. Figure [3] shows that the p ~ 3/2 limit is 
indeed realized if the MF is sufficiently broad and light- 
dominated. It corresponds to mass segregation stronger 
even than the saturation value pn = 5/4 predicted by 
AH09 in the limit of a strongly bimodal, light-dominated 
MF where light stars are expected to assume pl = 1/4. 

A variant of Eq. (f20|) with Q eliminated is 
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Adopting a typical specific energy, x^^-^^^^ with < 1, one 
can solve this equation for p{m) with arbitrary ^(m), an 
approximate procedure far simpler than solving Eq. ([9]). 
However, Eq. (|2ip typically becomes unstable when the 
underlying power-law assumption / ~ fails. 

The linear scaling p{m) oc m remains approximately 
valid even when the flow is only marginally negligible. 
Using this Ansatz in Eq. (pij). p{m) becomes a root of 
the quadratic equation 

4p^(m^) — p{7 — 4:p)m{m ) + {3/2 -p)m^ = , (22) 

with averages weighted by m£,x^^^. For light-dominated 
MFs and massive stars with m ^ {{'m), (m^)^/^}, the 
last term dominates and p{m) peaks at 3/2. For a 
strongly peaked MF where (m^) ~ {my, Eq. (|^^ yields 
two solutions: the Q = limit p = m/4(m), and a steep 
solution p = (3/2)(l + {m)/m)~^. 

Eq. (|22p suggests the following, numerically supported 
picture. For MFs peaked at relatively massive stars, the 
DF is approximately given by the Q = Q result. For 
more light-dominated MFs, Mq is shifted towards lower 
masses so the DF becomes gradually steeper at any given 
mass. If the MF is sufficiently extended (large C) and 
light-dominated {e.g., small a), each stellar species of 
mass m ^ eventually achieves maximal steepness 
p{m) ~ 3/2. The transition between the Q = and the 
saturation regimes involves in general an intermediate 
range of parameters in which the DF of this species is no 
longer a power law in x. We may roughly identify this 



region with the absence of real solutions to Eq. ([22]) . 

5{m) = {m/{m) - 5)^ - 24:{{m'^) / {m)^ - 1) < , (23) 
illustrated as areas delimited by solid curves in Figure 
[3l After reaching maximal steepness f{m) continues to 

evolve, but its maximal value remains ~ f{x — l)xmax- 

4. CONCLUSIONS 

We have analyzed the distribution of an arbitrary, con- 
tinuous MF ^ of stars around a MBH. For equal mass 
stars, we derive a simple approximation for the DF /, see 
Eq. ((8]) and Figure [T] For MFs not strongly dominated 
by light stars, where the mass flow can be neglected, we 
generalize the BW77 linear scaling p oc m to continu- 
ous MFs. We then derive approximate solutions for p 
(including its normalization, Eq. (fTS)) ') and / (Eq. (fT6|) ) 
and confirm them numerically; see Figures [2] and [3l 

Our results provide a simple yet accurate alternative to 
solving the full integro-differential equation ^ . Eq. (|15p 
reproduces and generalizes previous Fokker-Planck cal- 
culations such as BW77, which were limited to a narrow 
range of parameters. As an illustration, consider SgA* 
with Xmax — 10^ and a model MF with main sequence 
stars, white dwarfs, neutron stars and black holes, of 
masses M/Mq = 1,0.6, 1.4, 10, and relative abundances 
1 : 10-1 : 10-2 : IQ-^. Eq. ^ then gives Mq = A.SMq 
so ph — 0.52 for black holes and p < 0.08 for the 
ot her species, in agreement with the numerical results 
of iHopman and Alexande^ ()2006f ) and AH09. 

The DF becomes gradually steeper with decreasing 
Mq, as long as 5 in Eq. remains positive. For very 
light-dominated MFs, the power law assumption / ~ ^x^ 
eventually fails at high masses, roughly where S{m) < 
(contour-enclosed regions in Figure [3]). As the MF peak 
shifts to even lower masses, S{m) > and the power- 
law behavior are eventually restored at the high mass 
end and the DF remains very steep, peaking at p = 3/2. 
This behavior is illustrated in Figure [3] 

Power law DFs f (x x^ correspond to n{r) cx r^^/^^P. 
As long as the MF is not very light-dominated, stars of 
mass Mq have po = 1/4 (up to a p correction) and an 
equal-mass like cusp no cx r"^/"*, while other stars have 
nm/no cx r-(™-A^o)/4Mo^ jf ^.j^g j^p jg ^^^.^g^^ g^j^^ heavy- 
dominated, light stars tend towards an energy indepen- 
dent distribution n cx r^^/^ while heavy stars approach 
n cx r"^/*. In the opposite, maximally steep limit, heavy 
stars develop a cusp as steep as n cx , harder than any 
reported previously, and the number of stars in the cusp 
depends (logarithmically) on the disruption cutoff. 
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